QC using RSEM data

counts <- read_csv("data/genes_counts.csv")
tpm <- read_csv("data/genes_tpm.csv")
treatments <- read_csv("data/treatments.csv", col_types = "ccc")
ercc <- left_join(read_tsv("data/ercc.tsv"), tpm, by = "gene_id") 
descriptions <- read_tsv("data/descriptions.tsv", col_names = c("gene_id", "description"))
mix1 <-  ercc[,c("gene_id", "mix1", subset(treatments, spikein == "mix1")$sample)] %>% gather( library, observed, -c(mix1,gene_id)) %>% dplyr::rename(expected = mix1 )
mix2 <- ercc[,c("gene_id", "mix2", as.character(subset(treatments, spikein == "mix2")$sample))] %>% gather( library, observed, -mix2, -gene_id) %>% dplyr::rename(expected = mix2 )
ercc <- rbind(mix1,mix2)
ggplot(ercc, aes(observed+1, expected+1, color=as.factor(library))) + geom_point() + scale_y_continuous(trans='log2')+scale_x_continuous(trans='log2')+ geom_abline()+guides(color=F)

cutoff <- seq(0,5,by=0.1)
r <- c()
for (i in cutoff) {
  keepGenes <- ercc %>% group_by(gene_id, library) %>% summarise(all_mean = mean(observed)) %>% filter(all_mean > i) %>% .$gene_id %>% unique
  r <- c(r, summary(lm(log2(observed+1)~log2(expected+1), data=subset(ercc, gene_id %in% keepGenes)))$r.squared)
}
plot(cutoff,r)

The signal to noise is basically fine even at the lowest levels, so we do not do any abundance-based filtering.

Weighted gene co-expression analysis

We are interested in understanding the netwoks of genes co-expressed with Dsx in honey bees, and to checke whether this network is associated with RNAi treatment.

expDat <- as.data.frame(t(tpm %>% filter(! grepl("ERCC", gene_id)) %>% dplyr::select(-gene_id)))
colnames(expDat) <- tpm %>% filter(! grepl("ERCC", gene_id)) %>% .$gene_id
expDat <- expDat[as.character(treatments$sample),] # re-arrange to match treatments vector
pickSoftThreshold(expDat, RsquaredCut = 0.8, networkType = "unsigned")$powerEstimate

##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1    0.645  8.410          0.956  3900.0    3990.0   6160
## 2      2    0.163  0.263          0.602  1770.0    1790.0   3600
## 3      3    0.175 -0.348          0.817   985.0     970.0   2380
## 4      4    0.469 -0.672          0.845   620.0     579.0   1730
## 5      5    0.598 -0.824          0.866   423.0     370.0   1330
## 6      6    0.687 -0.902          0.893   305.0     249.0   1050
## 7      7    0.733 -0.947          0.896   230.0     173.0    857
## 8      8    0.762 -0.995          0.894   179.0     124.0    721
## 9      9    0.769 -1.050          0.876   142.0      91.3    623
## 10    10    0.776 -1.090          0.866   116.0      68.3    545
## 11    12    0.889 -1.090          0.945    80.7      40.2    431
## 12    14    0.898 -1.150          0.952    59.1      24.9    362
## 13    16    0.896 -1.230          0.954    44.8      16.1    318
## 14    18    0.894 -1.290          0.955    35.1      10.7    283
## 15    20    0.910 -1.300          0.971    28.1       7.4    255

## [1] 12

SoftPower = 12

We’ll pick soft power 12, which is over the 0.8 threshold.

net <- blockwiseModules(expDat, power = SoftPower,
TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs = FALSE,
verbose = 3, maxBlockSize = 15000)
saveRDS(net, "data/net.rds")

MEs0 <- moduleEigengenes(expDat, net$colors)$eigengenes
MEs <- orderMEs(MEs0)
modNames <- substring(names(MEs), 3)
moduleTraitCor <- cor(MEs, as.numeric(gsub("GFP",0,gsub("DSX",1,treatments$treatment))), method = "spearman");
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, 20)
(DsxModule <- net$colors[which(colnames(expDat) == "Dsx")])

## [1] 6

Alldegrees <- intramodularConnectivity.fromExpr(expDat, net$colors, power = SoftPower)

##  softConnectivity: FYI: connecitivty of genes with less than 7 valid samples will be returned as NA.
##  ..calculating connectivities..

rownames(Alldegrees) <- colnames(expDat)
geneModuleMembership <- as.data.frame(cor(expDat, MEs, use = "p"));

## Warning in cor(expDat, MEs, use = "p"): the standard deviation is zero

Alldegrees["Dsx",]

##       kTotal  kWithin     kOut     kDiff
## Dsx 76.54159 32.44871 44.09288 -11.64418

DsxModuleGenes <- colnames(expDat)[net$colors == DsxModule]
connectivityDsxModule <- data.frame(kWithin = Alldegrees[DsxModuleGenes,"kWithin"], gene_id = rownames(Alldegrees[DsxModuleGenes,]), moduleMembership = geneModuleMembership[DsxModuleGenes, paste0("ME",DsxModule)]) %>% left_join(de.tgw$table, by = "gene_id")

## Warning: Column `gene_id` joining factor and character vector, coercing
## into character vector

dsxModuleGenes <- descriptions %>% filter(gene_id %in% colnames(expDat)[net$colors == DsxModule] ) %>% left_join(de.tgw$table, by = "gene_id") %>% left_join(connectivityDsxModule %>% dplyr::select(gene_id, kWithin), by="gene_id") %>% dplyr::select(gene = gene_id, description = description.x, logFC, kWithin, logCPM) %>% arrange(desc(logFC))  
highligh <- which(dsxModuleGenes$gene %in% c("Vg", "Dsx", "Dnmt3"))